
*Welfare analysis simulation based on matlab results and the existing two policies (all combinations)

clear all

*read in the data: 

use "$root/Data/Produced/EV_welfare04.dta", clear

order price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4, after(choice)


*For belief error to be properly calculated var_costs need to be zero, where they are missing: 

replace var_cost=0 if missing(var_cost)

sort id car_option

*Set up the matrices and vectors needed for the computation:
mata

ndraws=150

b_f= (0.021183645,-2.593121925,0.154611204,-0.321820168,0.037136405,0.167413791,0.378321531,0.614503527,0.858999103,-0.005445116,-0.010654642,0.136554812,-0.083769393,-0.022510516,0.007285076,0.06246904,1.476449254,2.260850432,0.424973893,0.853075696,0.6579357,1.4180835,0.027265514,0.980066634,0.090419384,0.171157248,-0.142510175,-0.661616614,-0.177476937,0.595728981,0.214934965,0.378106038,-0.080343758,0.394597341,0.734410779,0.078025015,1.019999496,0.291640034,-0.620171333,-0.126470184,-0.114005609,-0.169462474,0.335901857,-0.232550483,0.070071269,0.05306741,0.003621919,0.006019563,0.024778481)'
/*
** Define the estimated vectors of random coefficients:
*/
beta_f=J(1,ndraws,b_f)
rseed(1234)
beta_p = rnormal(1,ndraws,-0.080965469,0.013364931)
beta_he = rnormal(1,ndraws,0.263640717,0.000661958)
beta_we = rnormal(1,ndraws,0.000192522,0.000045627)
beta_var = rnormal(1,ndraws,-0.040212482,0.001878464)


beta= (beta_p \ beta_he \ beta_we \ beta_var \ beta_f)

mata drop beta_he beta_we b_f


wealth_group=.
st_view(wealth_group,.,"wealth_quart2 wealth_quart3 wealth_quart4")


/*
Define the price coefficients:
*/
b_p=(0.003621919,0.006019563,0.024778481)'


b_w=wealth_group*b_p

mata drop wealth_group b_p



b_price=J(rows(b_w),1,beta_p) + J(1,ndraws,b_w)

gamma=beta_var:/b_price

mata drop beta_p b_w beta_var

end

timer on 1
mata
X=.
st_view(X,.,"price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_wealth2 price_wealth3 price_wealth4")


l_fit=exp(X*beta)

var_cost=.
st_view(var_cost,.,"var_cost")



/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id")
V=panelsetup(id, 1)

    ncol   = ndraws
    bigtot = J(rows(id), ncol, .)
    for (i = 1; i <= rows(V); i++) {
        X1 = panelsubmatrix(l_fit, i, V)
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1))
}

mata drop X1 ncol V

_editmissing(l_fit,0)

prob_all=l_fit:/bigtot
mata drop l_fit X 

prob=mean(prob_all')

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

mata drop prob_all var_cost be

st_addvar("double", "prob")
st_store(.,"prob",prob')

st_addvar("double", "be_init")
st_store(.,"be_init",be_m')


CS_init=mean((((-1):/b_price):*log(bigtot))')


mata drop bigtot

st_addvar("double", "CS_init")
st_store(.,"CS_init",CS_init')


mata drop prob CS_init be_m
end
timer off 1
timer list 1

*Generate the necessary frames to calculate the results: 
frame create results_pol subsidy tax_sim co2_ann emb_co2 petrol_tax tax_total tax_pv veh_tax paid_subsidy mean_co2 prob_EV 
frame create results_CS subsidy tax_sim CS_change be_change CSexp_change weight_CS_d weight_be_d weight_CSexp_d weight_CS_d_sq weight_be_d_sq weight_CSexp_d_sq alt_km
frame create results_driving subsidy tax_sim KM_driven_average KM_driven_obs

** Create the loop to variate both the subsidy and the tax policy: 

gen elast=-0.3

forvalues x=20/29 {
gen car_tax_alt=car_tax
replace car_tax_alt=car_tax/0.4 if drivetype==3
replace car_tax_alt=car_tax/0.6 if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax/0.8 if drivetype!=3 & env_cat==2	
	
	gen tax_sim=`x'
	
	merge m:1 tax_sim using "$root/Data/Produced/tax_params_ext.dta", 
	drop if _merge==2
	drop _merge
	
	
	
local r=0.06

gen pv_tax_alt=pv_tax

*for EV
replace pv_tax_alt=car_tax_alt*(1+EV_K)+(car_tax_alt*(1+EV_K))/((1+`r')^(1))+(car_tax_alt*(1+EV_K))/((1+`r')^(2))+(car_tax_alt*(1+EV_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype==3 

*for category A
replace pv_tax_alt=car_tax_alt*(1+A_K)+(car_tax_alt*(1+A_K))/((1+`r')^(1))+(car_tax_alt*(1+A_K))/((1+`r')^(2))+(car_tax_alt*(1+A_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

*for category B
replace pv_tax_alt=car_tax_alt*(1+B_K)+(car_tax_alt*(1+B_K))/((1+`r')^(1))+(car_tax_alt*(1+B_K))/((1+`r')^(2))+(car_tax_alt*(1+B_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*for category E
replace pv_tax_alt=car_tax_alt*(1+E_K)+(car_tax_alt*(1+E_K))/((1+`r')^(1))+(car_tax_alt*(1+E_K))/((1+`r')^(2))+(car_tax_alt*(1+E_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==5
*for category F
replace pv_tax_alt=car_tax_alt*(1+F_K)+(car_tax_alt*(1+F_K))/((1+`r')^(1))+(car_tax_alt*(1+F_K))/((1+`r')^(2))+(car_tax_alt*(1+F_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==6
*for category G
replace pv_tax_alt=car_tax_alt*(1+G_K)+(car_tax_alt*(1+G_K))/((1+`r')^(1))+(car_tax_alt*(1+G_K))/((1+`r')^(2))+(car_tax_alt*(1+G_K))/((1+`r')^(3))+(car_tax_alt)/((1+`r')^(4))+(car_tax_alt)/((1+`r')^(5))+(car_tax_alt)/((1+`r')^(6))+(car_tax_alt)/((1+`r')^(7)) +(car_tax_alt)/((1+`r')^(8))+(car_tax_alt)/((1+`r')^(9))+(car_tax_alt)/((1+`r')^(10))+(car_tax_alt)/((1+`r')^(11))+(car_tax_alt)/((1+`r')^(12))+(car_tax_alt)/((1+`r')^(13))+(car_tax_alt)/((1+`r')^(14)) if drivetype!=3 & env_cat==7

replace car_tax_alt=car_tax_alt*(1+G_K) if drivetype!=3 & env_cat==7
replace car_tax_alt=car_tax_alt*(1+F_K) if drivetype!=3 & env_cat==6
replace car_tax_alt=car_tax_alt*(1+E_K) if drivetype!=3 & env_cat==5
replace car_tax_alt=car_tax_alt*(1+B_K) if drivetype!=3 & env_cat==2
replace car_tax_alt=car_tax_alt*(1+A_K) if drivetype!=3 & env_cat==1
replace car_tax_alt=car_tax_alt*(1+EV_K) if drivetype==3 

gen ann_cost=drive_cost_ann+car_tax
gen ann_cost_alt=drive_cost_ann+car_tax_alt

gen alt_km=average_use*(1+elast*((ann_cost_alt-ann_cost)/ann_cost))

gen alt_cost_drive=(alt_km/100)*drive_cost
gen alt_drive_cost=alt_cost_drive+alt_cost_drive*((1-((1+`r')^(-14)))/`r')

gen alt_var_cost4=(alt_drive_cost+pv_tax_alt)/1000

drop alt_drive_cost alt_cost_drive ann_cost ann_cost_alt

replace alt_var_cost4=0 if missing(alt_var_cost4)

forvalues i=0 (200) 10000 {

timer on 1 
	
		gen subsidy=`i'
	
	gen price_alt_tot=prhnewprice
replace price_alt_tot=prhnewprice-subsidy if drivetype==3

gen price_alt=price_alt_tot/1000

forvalues a=2/4 {
	
	gen price_w_alt`a'=price_alt * wealth_quart`a'

}

drop price_alt_tot

sort id car_option

*Add the mata program here: 

mata {
X=. ;

st_view(X,.,"price_alt height weight alt_var_cost4 typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly resid price_w_alt2 price_w_alt3 price_w_alt4") ;

l_fit=exp(X*beta) ;

var_cost=.
st_view(var_cost,.,"alt_var_cost4")


/*
for (i=1; i<=ndraws; i++) {
st_addvar("double", "fit_" + strofreal(i))
st_store(.,"fit_" + strofreal(i),l_fit[,i])
}
*/
id=st_data(.,"id") ;
V=panelsetup(id, 1) ;

    ncol   = ndraws ;
    bigtot = J(rows(id), ncol, .) ;
    for (i = 1; i <= rows(V); i++) { ;
        X1 = panelsubmatrix(l_fit, i, V) ;
        bigtot[|V[i,1], 1 \ V[i, 2], ncol|] = J(rows(X1), 1, colsum(X1)) ;
} ;

/*
mata drop X1 ncol V
*/
X1=. ;
V=. ;
_editmissing(l_fit,0) ;

prob_all=l_fit:/bigtot ;
/*
mata drop l_fit X 
*/
l_fit=. ;
X=. ;

prob=mean(prob_all') ;

be=prob_all :* (J(1, cols(gamma), var_cost)-(gamma:*var_cost))

be_m=mean(be')

/*
mata drop prob_all
*/
prob_all=. ;
var_cost=. ;
be=. ;
st_addvar("double", "prob_alt") ;
st_store(.,"prob_alt",prob') ;

st_addvar("double", "be_alt") ;
st_store(.,"be_alt",be_m') ;



CS=mean((((-1):/b_price):*log(bigtot))');

/*
mata drop bigtot
*/
bigtot=.;
st_addvar("double", "CS_alt");
st_store(.,"CS_alt",CS');

/*
mata drop prob CS
*/
prob=.;
CS=.;
be_m=. ;
}


*Create a frame to calculate the consumer welfare change:
frame put CS_alt CS_init be_alt be_init id tax_sim subsidy TOTVERM prob_alt alt_km, into(cons)
frame put alt_km prob_alt choice id tax_sim subsidy, into(driving)

frame put prob_alt et_co2 alt_km et_verbrauch el_verbrauch car_tax_alt pv_tax_alt drivetype subsidy tax_sim car_option prhnewprice, into(probs)

frame change probs

gcollapse (mean) prob_alt et_co2 alt_km et_verbrauch el_verbrauch car_tax_alt pv_tax_alt prhnewprice (firstnm) drivetype subsidy tax_sim, by(car_option)

gen co2_ann=prob_alt*9230*et_co2*alt_km/1000

gen emb_co2_ann=0
replace emb_co2_ann=prob_alt*9230*alt_km/100*et_verbrauch*0.523 if drivetype==1
replace emb_co2_ann=prob_alt*9230*alt_km/100*et_verbrauch*0.445 if drivetype==2
replace emb_co2_ann=prob_alt*9230*alt_km/100*el_verbrauch*0.139 if drivetype==3
replace emb_co2_ann=prob_alt*9230*alt_km/100*el_verbrauch*0.139 + prob_alt*9230*alt_km/100*et_verbrauch*0.523 if drivetype==4


gen petrol_usage=prob_alt*9230*alt_km/100*et_verbrauch
gen petrol_tax=petrol_usage*0.7312
replace petrol_tax=petrol_usage*0.7587 if drivetype==2

gen tax_total=prob_alt*9230*car_tax_alt
gen tax_pv=prob_alt*9230*pv_tax_alt


gen veh_tax=prob_alt*9230*prhnewprice*0.04
replace veh_tax=0 if drivetype==3

gen paid_subsidy=prob_alt*9230*subsidy if drivetype==3
gen prob_EV=prob_alt if drivetype==3


gen mean_co2=prob_alt*et_co2


gcollapse (sum) co2_ann emb_co2_ann petrol_tax tax_total tax_pv veh_tax paid_subsidy prob_EV mean_co2 (firstnm) subsidy tax_sim


frame post results_pol (subsidy) (tax_sim) (co2_ann) (emb_co2_ann) (petrol_tax) (tax_total) (tax_pv) (veh_tax) (paid_subsidy) (mean_co2) (prob_EV)

frame change cons

*Collapse in case two options had the exact same predicted prob: 
gcollapse (firstnm) CS_alt CS_init subsidy tax_sim TOTVERM alt_km (sum) be_alt be_init , by(id)

*Generate the CS_change
gen CS_change=CS_alt-CS_init
gen be_change=be_alt-be_init
gen CSexp_change=CS_change-be_change

*Rescale so it roughly represents one year of registrations:
replace CS_change=9230/23074*CS_change
replace be_change=9230/23074*be_change
replace CSexp_change=9230/23074*CSexp_change

*Calculate the Welfare function of the social planner
gen weight_CS_d=(1/TOTVERM)*CS_change
gen weight_CSexp_d=(1/TOTVERM)*CSexp_change
gen weight_be_d=(1/TOTVERM)*be_change

gen weight_CS_d_sq=(1/(TOTVERM^2))*CS_change
gen weight_CSexp_d_sq=(1/(TOTVERM^2))*CSexp_change
gen weight_be_d_sq=(1/(TOTVERM^2))*be_change


gcollapse (firstnm) subsidy tax_sim (sum) CS_change be_change CSexp_change weight_CS_d weight_be_d weight_CSexp_d weight_CS_d_sq weight_be_d_sq weight_CSexp_d_sq  alt_km

frame post results_CS (subsidy) (tax_sim) (CS_change) (be_change) (CSexp_change) (weight_CS_d) (weight_be_d) (weight_CSexp_d) (weight_CS_d_sq) (weight_be_d_sq) (weight_CSexp_d_sq) (alt_km)

frame change driving
gen drive_obs=choice*alt_km
gen drive_prob=prob_alt*alt_km
gcollapse (sum) drive_obs drive_prob (firstnm) subsidy tax_sim


frame post results_driving (subsidy) (tax_sim) (drive_prob) (drive_obs)

frame change default	

frame drop cons
frame drop probs
frame drop driving

drop subsidy prob_alt CS_alt price_alt price_w_alt2 price_w_alt3 price_w_alt4 be_alt

timer off 1
timer list 1

}

drop tax_sim alt_var_cost4 car_tax_alt pv_tax_alt A_K B_K E_K F_K G_K EV_K alt_km

}


**After the calculations bring the result frames together to calculate the welfare outcomes: 
frame change results_pol

frlink 1:1 subsidy tax_sim, frame(results_CS)
frlink 1:1 subsidy tax_sim, frame(results_driving)
frget _all, from(results_CS)
frget _all, from(results_driving)

drop results_CS results_driving


** Save the datafile
save "$root/Results/estimation_results/welfare_max08_3y_elast_3.dta", replace


